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A  semiclassical  method  for  solving  the  quantum  Liouville  equation  in 
one-dimensional  phase-space  is  described*  The  development  is  based  on 
constructing  a  Gaussian  density  matrix  and  is  applicable  to  systems  in  pure 
and  in  mixed  states  having  nonlinear  interaction  potentials.  The  density 
matrix  is  constructed  using  a  set  of  dynamic  variables  whose  expectation 
values  are  considered  to  be  relevant  for  the  dynamics.  The  self-consistent 
equations  of  motion  are  then  derived  for  these  expectations  from  the  quantum 
Liouville  equation  using  a  projection  scheme-  The  solution  of  these  self- 
consistent  equations  provides  the  time  evolution  of  the  density  matrix.  The 
present  method  can  yield,  in  principle,  exact  values  for  these  expectations 
for  all  times.  A  model  calculation  is  carried  out  to  describe  the 
vibrational  motion  of  an  arbitrary  diatomic  molecule  on  an  anharmonic 
potential  surface.  However,  the  potentiality  of  this  method  lies  in 
describing  the  time  evolution  of  systems  in  mixed  states  and  hence  in 
describing  the  dynamics  of  molecular  processes  in  condensed  phases. 


I.  INTRODUCTION 


Recent  advances  in  the  experimental  study  of  the  various  molecular 
dynamical  processes  in  condensed  phase,  such  as  energy  transfer,  molecular 
dissociation  reactions,  spectral  line  shapes,  etc.,  require  theoretical 
models  for  the  quantitative  understanding  of  the  dynamical  processes  involved 
in  condensed  phases.  There  has  been  progress  in  studying  equilibrium 
properties  using  classical  [1],  semiclassical  [2],  fully  quantum  mechanical 
[3]  and  quantum  field  theoretic  methods  [4,5].  Methods  are  also  available 
for  treating  time -dependent  processes  within  the  classical  framework  [6]. 
However,  very  few  theoretical  developments  are  available  for  treating  time- 
dependent  processes  incorporating  quantum  effects.  These  are  the  quantum 
corrections  to  the  classically  computed  time-correlation  functions  [7],  the 
exp(S)  approach  of  Arponen  ani  co-workers  [4]  and  the  semiclassical  Gaussian 
wavepacket  dynamics  (GWD)  approach  developed  notably  by  Heller  [8].  The 
semiclassical  GWD  approach  describes  a  self-consistent  solution  of  the  time- 
dependent  Schrodinger  equation  and  thus  is  restricted  in  its  application  to 
systems  in  pure  states.  Extension  of  this  GWD  method  to  the  simulation  of 
time-dependent  properties  of  N-particle  systems  interacting  through  realistic 
pair  potentials  within  the  variational  and  nonvariational  framework  are  also 
available  in  the  literature  [9J.  Such  application  requires  tedious  thermal 


averaging,  which  arises  from  the  fact  that  we  have  no  knowledge  about  initial 
conditions  of  the  N-particle  system. 

Our  objective  is  to  develop  a  similar  GWD  approach  which  as  such  is 
applicable  to  systems  in  pure  and  in  mixed  states.  That  is,  when  treating 
systems  in  mixed  states,  we  do  not  need  to  perform  tedious  thermal  averaging. 
Our  development  satisfies  the  maximum  entropy  principle  (10]  when  treating 
equilibrium  or  nonequilibrium  systems.  However,  we  no  not  make  the 
assumption  that  the  exact  nonequilibrium  statistical  density  matrix  is 
approximately  equal  to  the  local  equilibrium  one  [11-13].  For  an  N-particle 
statistical  system  it  is  practically  impossible  to  construct  a  density  matrix 
which  contain  all  information  about  the  system.  However,  with  the  advent  of 
projection  operator  techniques  [11,14],  it  has  been  possible  to  construct 
density  matrices  which  contain  information  sufficient  for  the  calculation  of 
various  physical  quantities  of  interest. 

In  this  paper  we  are  interested  in  a  reduced  description  of  the  exact  N- 

particle  system,  which  is  the  time  evolution  of  the  N  single-particle  density 

matrices  in  a  mixed  state.  We  define  our  reduced  density  matrix, 

p  (X,X',t),  as  a  product  of  N  single-particle  density  functions 
re 

N 

Pre(X,X';t)  -  ♦jUj.xjst)  ,  (1) 

where  X  is  a  vector  with  N  coordinate  components  x^...x^.  The  time  evolution 
of  these  density  functions,  are  then  obtainedAf rom  the  quantum  Liouville 
equation  using  a  projection  operator  scheme  [11,14].  We  define  each  single¬ 
particle  density  function  $.(x.,x!;t)  from  the  perspective  of  nonequilibrium 
statistical  mechanics  as  [lO^ll^lAi 

♦  »  <x  |i  (t>|«’> 

*  <Xj|expl£  Aja(t)  A^aJ |xj>  ,  (2) 

a-0 

which  contains  all  information  about  the  single  particle  system.  The 
A.  (t)'s  are  Lagrange  multipliers,  and  the  A.  *s  are  the  dynamical  variables. 
Since  we  are  not  interested  in  all  the  information  contained  in  the  ♦  's,  *e 
construct  our  +j,s  vith  respect  to  the  six  dynamical  quantities  ^ 

Ajo  "I*  *J1  "  *J2  "  *J*  *J3  "  *J ’ 

*J4  ■  4  ‘js-Vj'Vj  •  (J> 

where  is  the  momentum  associated  with  the  k-th  particle  and  the  hat 

designates  an  operator.  As  we  shall  see  later,  the  choice  of  these  dynamical 
quantities  allows  us  to  describe  the  time  evolution  of  each  single-particle 
density  function  incorporating  quantum  fluctuations.  The  time  evolution  of 
the  expectations  of  these  dynamical  quantities,  <A.  >»  are  then  obtained  from 
the  quantum  Liouville  equation  using  the  projectionoperator  scheme  [11,14]. 
The  choice  of  the  single-particle  density  operator  as  given  by  Eq.  (2)  is  by 
no  means  unique.  Our  choice  is  motivated  by  the  physical  consideration  which 
is  the  maximum  entropy  principle  [10,11]. 

We  confine  our  development  to  one-dimensional  phase  space.  In  the  next 
section  we  derive  the  equations  of  motion  for  the  expectations,  <A.  >,  in 
closed  form  and  construct  the  corresponding  density  function  forJ  mixed 
states.  In  Sec.  Ill  we  show  that  under  certain  conditions  the  density 


function  for  mixed  states  reduces  to  the  density  function  for  pure  states. 
To  describe  the  time  evolution  of  the  pure  state  density  function,  we  then 
derive  the  equations  of  motion  for  the  corresponding  dynamical  quantities. 
In  Sec.  IV  we  show  that  our  maximum  entropy-based  density  function  can  also 
describe  the  time  evolution  of  a  harmonic  system  in  thermal  equilibrium  1 15]. 
In  Sec.  V  we  solve  the  equations  of  motion  for  the  pure  state  to  describe  the 
vibrational  motion  of  an  arbitrary  diatomic  molecule  on  a  Morse  potential 
surface.  Vie  then  compare  our  results  with  those  obtained  using  the  existing 
Gaussian  wavepacket  dynamics  method  [8,9],  and  a  discussion  is  provided  in 
Sec.  VI. 

II.  CONSTRUCTION  OF  THE  DENSITY  FUNCTION  AND  DERIVATION  OF  THE  EQUATIONS  OF 
MOTIONS  FOR  SYSTEMS  IN  MIXED  STATES 

We  characterize  our  N-particle  system  by  a  Hamiltonian 
N  2 

»-£  ST  +  <4) 

k*l  k 

and  a  density  matrix  p(X,X';t)  which  satisfies  the  quantum  Liouville  eqution 

-  -iLp  5  -(i/K)[H,p]  ,  (5) 

where  m,  is  the  mass  of  the  k-th  particle  and  V  is  the  interaction  potential. 
Since  we  are  interested  only  in  the  time  evolution  of  the  N  single-particle 
density  functions,  we  partition  our  total  density  matrix  as 

p(t)  -  Pr<(t)  +  Pir(t)  .  (6) 

where  p  (t)  is  a  reduced  description  of  the  N-particle  interacting  system 
and  isrCrepresented  by  a  product  of  N  single-particle  density  functions  as 
described  in  Eq.  (1).  p.  (t)  represents  the  irrelevant  degrees  of  freedom, 
since  it  does  not  contain  any  dynamical  degrees  of  freedom  of  any  single 
particle  in  the  coupled  N-particle  system,  but  rather  the  correlations 
between  single  particle  systems  produced  by  their  interaction. 

We  associate  entropy  S  with  our  system  by  using  the  relation  [10] 

S  •  -  k  Trpr#(t)  lnpre(t)  (7) 

where  k  is  Boltzmann's  constant.  We  maximize  entropy  subject  to  the 

constraints 

Trp  (t)  ■  1  (8a) 

and 

«.  (t)  »  <A .  (t)>  -  TrA.  p  (t)  *  TrA.  p(t)  ,  (8b) 

Ja  jo  ja  re  ja 

where  the  Aj^'s  are  the  6N  dynamical  variables  of  interest  to  us. 

We  now  derive  explicit  expressions  for  the  time  evolution  of  the 
expectations,  a.  (t),  using  the  time-dependent  projection  operator  scheme 
[14]  followed  in^?onstructing  the  maximum  entropy  distribution  of  the  reduced 
density  operator  p  (t)  in  one -dimens tonal  phase  space.  From  now  on  we  shall 
refer  to  these  equifions  of  motion  as  reduced  equations  of  motion  since  they 
describe  the  time  evolution  of  the  reduced  density  operator  Pre(t).  We  shall 


use  the  projection  operator  technique  in  Liouville  space  [11,14].  In  this 
space  ft  and  p(t)  can  be  written  as  |H>>  and  p(t)>>.  In  this  notation  Eq. 
(8b)  becomes 


AJo(t)  *  <<pre(t) I Ajo>>  5  «P<t>lAJa»  • 


For  each  degree  of  freedom  j,  we  now  define  a  6  *  6  matrix  with  elements 

oJpU)  -  «Ajal pre(t)Ajp»  5  TrtAjapre^Ajp^ 

a,p  -  0,1, ...5  (9) 

and  the  Liouville  space  projection  operator 
N  5 

«« - 1  l  ip«(t)*Ja»[DJ(t)i;;  «*J8i  ao 

j-1  a, p-o 

having  the  following  properties: 

a)  P(t)  Pit')  -  P(t);  t  >  t' 

PZ(t)  -  P(t)  (11 


b)  P(t)|p(t)»  »  |pre(t)» 

c)  «Aja|P(t)p(t)»  -  «AJa| preCt>» 

where  p(t)  -  and  p  (t)  * 

at  re 

d)  P(t)  |p(t)Ajp»  -  lpre^t^Ajg>>  • 


dpre(t) 


(11a) 

(lib) 


(11c) 

(lid) 


It  has  been  shown  in  a  separate  conmunication  [16]  that  the  properties  (11) 
can  easily  be  derived  using  the  definitions  (9)  and  (10).  P(t),  therefore, 

is  the  projection  operator,  since  it  reduces  the  exact  density  matrix  p(t) 
to  the  simpler  distribution  pf#(t). 

Let  us  now  assume  that  at  some  time  t  -  t' 

P(f)  -  Pr.(t')  .  (12) 

Using  this  assumption  and  introducing  the  complementary  projection 


Q(t)  -  1  -  P(t) 


as  shown  in  Ref.  16,  we  can  write  the  exact  reduced  equations  of  motion 

(REM)  for  the  e,  (t)'s  from  the  quantum  Liouville  equation  (S)  as 
J® 

•jo(t)  "  '  i<<AjJLlpre(t)>>  *  l  Mig(t»t')*jp(t)  ’  (U*> 

where  we  have  introduced  the  6  *  6  matrices 

w^(t,f)  -  “i«A.  |LQ(t)U(t,t*  )  lpr#(t'  )Ajg»  *  (14b) 

R^(t.f)  -  «Aja|U(t,t')|pr%(f)Aje»  .  (14c) 

and 


(14d) 


Here  U(t,t')  is  the  time  evolution  operator 
U(t,t')  -  exp[-iL(t-t')J  . 

Equation  (14d)  can  be  recast  in  matrix  notation: 
MJ(t,t')  *  WJ(t,f)[RJ(t,t')]'1  . 


( 14e) 


( 14f  ) 


Equation  (14)  describe  the  time  evolution  of  the  5N  dynamical  quantities  a. 
(j  *  1,...N;  a  ■  1,2,.. .5)  and  are  exact.  There  are  5N  nonlinear  coupl^S 

differential  equations  for  5N  unknown  a  (t).  In  these  equations  the  time 
derivative  of  a,  at  time  t  depends  onfall  a  at  the  same  time.  Note  that 
ve  assume  A  Q  toJbe  the  unit  operator,  anaP  normalization  requires  its 
expectation*1  value  to  be  independent  of  time,  a  *  1.  An  alternative 
derivation  of  Eqs.  (14)  is  also  possible  [14,16],  where  the  time  derivative 
of  a  at  time  t  depends  on  all  a.g  at  previous  time  t1  <  s  <  t,  and  the 
resulting  equations  are 


•Ja(t>  ‘  •1<<»jalt-l“re<t)» 

-  J  ds  <<Aja|LH(t ,s)Q(s)L| p^c(t)» 


where 


H(t,s)  «  e 


xp[-i  f 

J  n 


dt  Q(t)L] 


(15a) 


(15b) 


is  a  time-ordered  exponential.  Now  if  we  assume  that  condition  (12)  holds 
for  all  times,  then  Q(t)p(t)  *  0,  and  we  are  left  with  the  first  term  on  the 
right-hand  side  of  both  Eqs.  (14a)  and  (15a),  which  represents  a  mean  field 
time  evolution  of  the  N-particle  system.  The  second  terms  are  the 
correlation  terms  and  arise  from  the  fact  that  p(t)  *  p  (t)  for  all  times. 
If  we  retain  up  to  a  given  order  in  the  correlation  termsCin  Eqs.  (14a)  and 
(15a),  then  they  yield  different  approximations.  However,  in  this  paper  we 
are  interested  only  in  the  mean  field  time  evolution  of  the  N-particle 
system,  where  the  time  evoution  of  the  expectations  of  the  dynamical 
quantities,  A^,  are  given  by 

iJa(t)  -  i/|l  Tr{Aja(H,pre(t)]}  •  (16) 

For  our  convenience,  however,  we  evaluate  explicitly  the  time  evolution  of 
the  dynamical  quantities 


°J1  "  <XJ>*  °J2  "  <PJ> 

°J3  "  <*J>  *  <XJ>2,  °J4  "  <p2>  ‘  <pj>2 


J5 

given  by  [16] 


+  Pj*j>  -  2<Xj><Pj>]  , 


(17a, b) 
(17c, d) 
(17e) 


°jl  •  Vj  • 


(18a) 


(18b) 

(18c) 


°J2 

°j3 

*J5 


where 


-<VJ(X)>  , 

0j5/mj  ’ 

-<Vj'(X)>Oj5  , 

2^T  '  <V','(X)>a .-}  . 

m.  J  j  J 


V'  (X)  ■ 

VjW  3Xj  * 

V"(X)  -  ~  , 

J  8«j 


<Vj(X)>  -  |  dX  Vj(X)pre(X,X;t)  , 
N 

Pr@(X,X;t)  -  n  (>j(Xj,Xj}t) 


♦j(Xj,Xj}t) 


77 


i  f  (*r#u)  . 

*rr exp[_  2^ — ] 


J3 


J3 


(18d) 

(18e) 

(19a) 

(19b) 

(19c) 

(19d) 

(19e) 


<V',(X)>  is  expressed  in  a  similar  way  to  Eq.  (19c)  by  replacing  V'(X)  with 
V"fX).  Equations  (18)  are  the  time-dependent  self-consistent  field^  (TDSCF) 
equations.  For  each  particle  j  we  obtain  a  closed  set  of  five  equations, 
which  show  correct  self-consistent  behavior  in  any  potential  and  are  coupled 
to  each  other.  The  first  two  equations  in  (18)  express  Ehrenfest's  Theorem 
[17],  and  the  third  and  fourth  give  a  measure  of  the  uncertainty  in  position 
and  momentum  measurement  in  the  system.  The  fifth  equation  appears  only 
when  we  are  treating  systems  in  a  mixed  state.  For  systems  in  a  pure-state. 


"  (4°j3°J4  ’  *  (20) 

From  Eq.  (19)  we  find  that  for  a  successful  application  of  the  REM,  the 
choice  of  the  form  of  ♦.  is  crucial.  Our  particular  choice,  as  described  by 
Eqs.  (2)  and  (19e),  is  by  no  means  unique.  We  are  motivated  by  the  physical 
consideration  which  is  the  maximum  entropy  principle  [10].  Such  choice  for 
#.  connects  the  present  semiclassical  procedure  with  the  more  general 
problem  of  the  derivation  of  REM  in  nonequilibrium  statistical  mechanics 
[18,20]. 


In  the  following  we  show  that  the  choice  of  the  dynamical  quantities  as 
described  by  Eq.  (3)  produces  a  Gaussian  distribution  for  each  4,  in  one¬ 
dimensional  phase  space.  We  now  derive  explicitly  the  phase-space 
representation  (q,p)  for  one  degree  of  freedom.  The  proof  holds,  however, 
for  any  N  since  we  represent  our  reduced  density  function  p  (X,X*;t)  by  a 
product  of  N  single-  particle  density  functions  $  (x.,x*;t)  (£$.  (1)).  We 
therefore,  from  now  on,  choose  to  drop  the  subscript  j  and  replace  Eq.  (1) 
by 

o(x,x',t)  £  <x|o(t)|x’>  (21a) 


with 


(21b) 


5 

o(t)  =  exp[)  A  (t)A  ) 
a=0 

As  shown  in  Ref.  16,  the  Wigner  representation  [19]  of  the  density  operator, 
o(t),  may  be  written  in  the  form 

°w(q»p;t)  *  £  lap-y2/4]^exp[(62&+(>2a-Y6i*)/(4a3-Y2)] 

2  2 

*  exp[aq  +f5p  +ypq+6q+^p] 


with 


II 


dqdp  Ow(q,p,t)  »  1 


(22) 


and  the  corresponding  coordinate  representation  is  obtained  from  the 
transformation 


c(q+s,q-s;t)  *  dp  ow(q,p;t)exp[2ips/tf J 

Using  the  substitutions 

q  *  (x+x')/2,  s  -  (x-x')/2 


(23) 


(24) 


in  Eq.  (23),  we  obtain 

o(x,x';t)  »  ^"'r2~4ot^^^exPf ( 2p6-Y(>)2/4p(4<xp-Yr2 ) ] 

2 

*  exp[^(a-Jg)(x+x' )2  +  l(<S-2p)(x+x' ) 

+  ^(X'X')2  "  ^(x2'x’2)  '  2gK(x'x,)1  ’  (25) 

with  a,  g,  6  and  4  being  real  parameters,  which  may  be  expressed  in  terms 
of  Equations  (22)  and  (25)  accomplish  our  goal  of  expressing  the 

maximum  entropy  distribution  [Eq.  (21b)]  in  phase  space  (q,p)  and  in  the 
coordinate  representation  (x,x*).  However,  to  obtain  the  TDSCF  set  of 
equations  (18),  we  have  used  a  different  form  of  representation  of  these 
distribution  functions,  which  were  obtained  by  expressing  Eqs.  (22)  and  (25) 


terms  of  the  expectations  of  the  dynamical  quantities 
.  They  are  related  to  the  parameters  by 

described  by  Eq. 

o1(t)  -  <x(t)>  »  (236-y$)(y2-4cx$)  1 

(26a) 

02(f)  *  <p(t)>  *  (2a^-y6)(Y2*Aap)  * 

(26b) 

o^(t)  ■  <x2(t)>  -  <x(t)>2  ■  2p(Y2~^a^)  * 

(26c) 

o4(t)  ■  <p2(t)>  -  <p(t)>2  ■  2a(Y2-4a&)  1 

(26d) 

Oj(t)  *  <xp+px>  -  2<xXp>  *  -2y(y2*^ci3)  * 

(26e) 

a  *  C<V  P  *  ca3*  y  9  “CO5 

6  *  c^°2°6  *  2o1°4^#  *  “  c^°i°5  ‘  2o2°3^  * 


or 


where 


c  =  i(Y  -4ap)  = 


oc  -40-0, 
5  3  4 


Expressing  the  phase-space  density  function  o  (q,p,t)  of  Eq.  (22)  in  terms 
of  the  Oj(t)'s,  we  have 


ou(q,p;t)  =  - 5“T  exp{- 

ir(4 


; — ^  lo4(<*'°i)2  +  °3{p'°2)2 
4o3V°5 


-  o5(q-o1)(p-o2)l}  , 

and  the  corresponding  coordinate  representation  (Eq.  (25))  becomes 


o(x,x' ;t) 


exp(-o12/2o3) 


;xp[-  ^-(x+x1)2  +  —•  (x+x')  +  — =r — (x-x')2 


"  ®°3  '  2o3  '  4K2o3c 

+  4^*  (x2-x'2>  *  (<?1o5-2o2o3)(x-x1)]  .  (29) 

Thus,  the  particular  choice  of  the  dynamical  quantities,  as  depicted  in  Eq. 
(3),  generates  a  Gaussian  form  for  the  representation  of  the  corresponding 
single-particle  distribution  functions.  In  the  following  section  we  shall 
show  that  condition  (20)  reduces  these  mixed-state  density  functions  [Eqs. 
(28)  and  (29)]  to  that  of  pure  states.  We  shall  also  derive  the  REM  for 
pure  states. 

III.  REDUCED  EQUATIONS  OF  MOTION  (REM)  FOR  THE  PURE-STATE  DENSITY  FUNCTION 

Using  the  same  projection  scheme  as  in  Sec.  II,  a  self-consistent 
description  for  the  time  evolution  of  the  pure-state  density  function  can 
also  be  obtained.  Following  Heller  [8],  we  define  the  reduced  density 
function  for  pure  states  in  the  coordinate  representation  as 

oH(x,x';t)  *  exp ( - ^  {(x-xt)2  +  (x'-xt)2} 

+  g  a2((x-xt)2  -  (x'-xt)2}  +  ^  Pt(x-x' )]  ,  (30) 

where  the  parameters  o. ,  a»,  x  and  p  are  related  to  the  o  (t)'s  as 
follows:  1  L  1  t  i 


Oj(t)  -  xfc,  o2(t)  -  pt  , 

°3(t)  *  4^’  °4(t)  * 

o5(t)  -  ^  -  (4o3(t)oA(t)-J(2]1 


°  "  lai  +  a2 


(31a, b) 


(31c, d) 


Therefore,  oH(x,x*;t)  -  75==  exp(-  ^-^expt-  ^-(x^+x'2)  +  ^-(x+x1) 


(31e) 


i 


(32) 


io 


4fto 


5  ,  2  , 2  x 

(x  “ x 1  )  * 


20  (o1°5'2o2<J3)(x'x,)1 


and  the  corresponding  phase-space  representation  is  given  by 

°HW(q*p,t)  *  O  6Xpl' 

+  o3(p-o2)2  -  Oj(q-o1)(p-o2)}]  ,  (33) 

where  is  given  by  (31e),  which  is  the  same  as  condition  (20).  These 
pure-state  density  functions  can  also  be  obtained  directly  from  the  mixed- 
state  density  functions  (Eqs.  (28)  and  (29))  using  condition  (20). 


We  now  assume  that  the  time  evolution  of  the  pure-state  system  is 
described  by  the  approximate  density  functions,  (32)  and  (33),  for  all 
times.  This  assumption  then  allows  us  to  construct  the  SCF  set  of  equations 
for  the  expectations  of  the  corresponding  dynamical  quantities  using  Eq, 
(16).  They  are 

0-(t) 

°l(t)  =  m 

(34a) 

o2(t)  *  -<V’(x)> 

(34b) 

c>3(t)  *  ~  (4o3(t)o4(t)-H2)^ 

(34c) 

o^(t)  =  -(4o3(t)o4(t)-tt2)^  <V"(x)>  . 

(34d) 

This  is  a  closed  set  of  four  equations  which  differ  from 
equations  in  the  mixed  case  [Eqs.  (18a)  through  (18d)]  due 
is  no  longer  an  independent  variable  [Eq.  ( 20 ) ] ,  and  for 
we  do  not  have  any  REM  for  in  the  pure  case. 

the  first  four 
to  the  fact  that 
the  same  reason 

The  REM  of  Eq.  (34)  are  obtained  from  Eq.  (16)  under  the  exact 
potential  of  the  system.  These  equations,  as  shown  below,  are  different 
from  those  obtained  by  solving  the  quantum  Liouville  equation  (5)  for  pure 
states  under  the  locally  quadratic  potential  approximation  [8).  They  are 

0,(t) 

°l(t)  *  ffl 

(35a) 

o2(t)  -  -V(*)|x.0i 

(35b) 

<J3(t)  »  (4o3(t)o^(t)-h2)^/m 

(35c) 

o4(t)  -  -(4o3<t)o4(t)-#2)V(x)|  . 

(35d) 

Here  the  first  two  equations  describe  the  classical  motion  for  a  system  in  a 
pure  state  and  are  not  coupled  with  the  other  two  equations  [(35c)  and 
(35d)],  which  describe  the  time  evolution  of  the  variances  (Eqs.  (17c), 

( 1 7d ) ) .  Therefore,  the  present  set  of  equations  describe  the  trajectory  of 

a  particle  whose  position  and  momentum  at  time  t  are  known  from  the  center 
of  the  wavepacket.  However,  the  trajectory  of  a  particle  is  described  by 
the  wavepacket  [Eqs,  (32)  and  (33)1  as  a  whole,  which  inevitably  has  certain 
spacial  extension. 


On  the  other  hand,  if  we  look  at  the  SCF  set  of  equations  (34),  we  find 
that  the  first  two  equations  (a  and  b)  are  coupled  with  the  other  two 
equations.  Again,  the  right-hand  side  of  Eq.  (34b)  is  equal  to  the  average 
of  the  force  over  the  whole  wavepacket  and  thus  differs  from  Eq.  (35b)  due 
to  the  fact 

<V' (x)>  *  V* (x) I  .  ( 36) 

In  Sec.  V  we  shall  analyze  the  relative  meris  of  these  approaches  by 
studying  the  vibrational  motion  of  an  arbitrary  diatomic  molecule  on  a  Morse 
potential  surface.  It  is  important  to  note  here  that  for  the  mixed  case, 
even  if  we  start  with  a  minimum  uncertainty  wavepacket  the  variances  and 
are  not  constants  of  motion  as  the  system  evolves,  which  is  evident  from 
Eqs.  (18).  In  deriving  expressions  (32)  and  (33)  we  assumed  that  the  N- 
particle  density  function  may  be  written  as  [ 20 J 

p(X,X* ; t )  -  Y(X,t)  ¥*(X\t)  ,  (37) 


where  ?(X,t)  is  the  exact  wave  function  of  the  N-particle  interacting 
system.  We  then  introduced  the  approximation 


y(x,t) 


N 

*  n 

j=i 


(38) 


where  the  .  (x  ,t)'s  are  the  single-particle  wave  functions  and  contain  all 
information^ about  the  single-particle  systems,  including  their  phase.  A 
reduced  description  of  these  single-particle  wave  functions  was  first 
introduced  by  Heller  [8],  which  in  terms  of  the  o^'s  may  be  written  as 


♦  (x,t)  -  (2ho3)  ^  exp{(^~  + 


io-  2 

Hjc-Oj)  ♦  “jr(x’°l)  + 


ii) 

r 


(39) 


where  for  notational  convenience  we  have  dropped  the  subscript  j.  is 
given  by  condition  (20),  and  y  is  the  phase-f actor .  The  density  function 
corresponding  to  this  wave  function  is  given  by  expression  (32),  which  we 
obtained  from  the  maximum  entropy  distribution  (29)  using  condition  (20). 
Therefore,  if  we  assume  this  Gaussian  wave  function  (39)  to  approximate  the 
exact  single  particle  wave  function  for  all  times,  the  time  evolution  of 
this  reduced  wave  function  under  the  exact  potential  of  the  system  can  be 
obtained  by  solving  the  SCF  set  of  REM  given  by  Eqs.  (34),  along  with  the 
equation  for  the  phase  factor 


y(t) 


which  is  obtained  from  the  Schrodinger  equation 

2 

<E>  -  +  V(x)>  *  iK<*|*> 


(AO) 


(41) 


The  quantities  Vq, 


<V>,  V, 


in  Eq. 


(40)  are  given  by 


Heller  first  evaluated  this  propogator  (39)  under  the  locally-quadrat ic 
potential  approximation  [8J. 

IV.  CANONICAL  DENSITY  FUNCTION  FOR  A  HARMONIC  SYSTEM 


In  this  section  we  show  that  the  TDSCF  set  of  equations  (18),  which 
describe  the  time  evolution  of  any  irreversible  process  under  the  exact 
potential  of  the  system  by  using  the  reduced  density  matrix  expressions  (28) 
and  (29),  can  be  used  to  describe  the  time  evolution  of  a  harmonic  system  in 
thermal  equilibrium  [15].  When  a  system  is  in  thermal  equilbrium,  we  have 
the  density  matrix  satisfying  maximum  entropy  principle  as  [10] 

oT(t)  *  exp(-pH)/Tr[exp(-pH)]  ,  (42) 

where  g  *  (kT)  *  and  H  is  the  Hamiltonian  of  the  system.  Under  the 
quadratic  potential  approximation,  where 

H(q,p)  *  p2/2m  +  irow2q2  ,  (43) 

a  Gaussian  form  of  representation  of  the  density  operator  (42)  can  be 
obtained,  which  in  the  coordinate  representation  (x,x‘)  is  given  by 

oT(x.,.;t)  -  1°“ 

•  exP<2»  slnM8E;1(l,2t>'2)':osh<el>“)~2l'<'l>  '  (44) 

Expectation  values  of  the  dynamical  quantities  (17),  with  respect  to  this 
density  matrix,  are 


ox  -  0  ,  o2  -  0 

°3  *  9^  coth(if5Jta»)  ,  o4  *  Jmwji  coth(ig)ta)  ,  *  0  ,  (45) 


where  for  convenience  we  have  dropped  the  j -subscript.  Now  expressing  the 
thermal  density  function  in  terms  of  the  o^(t)'s,  we  obtain 


oT(x,x' ;t) 


===  exp{-  t^-(x+x')2 - ^r(x-x')2}  , 

xo3  803  n2 


and  the  corresponding  phase-space  density  function  becomes 


°TV(q’PJt)  *  exp(-  2? 
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The  time  evolution  of  these  density  functions  are  found  by  solving  the  set 
of  equations  (18)  with  initial  conditions  given  by  Eq.  (45)  and  the 
interaction  potential  given  by  expression  (43).  This  is  because  Eqs.  (46) 
and  (47)  do  not  satisfy  condition  (37).  Our  development,  as  described  in 
Sec.  II,  however,  is  more  general  since  it  can  be  used  for  studying  the 
relaxation  of  a  system  to  thermal  equilibrium  with  a  thermal  bath  under  the 
exact  potential  of  the  system. 


V.  VIBRATIONAL  MOTION  OF  AN  ARBITRARY  DIATOMIC  MOLECULE  ON  A  MORSE 
POTENTIAL  SURFACE 


In  this  section  we  solve  the  TDSCF  set  of  RQi  (34)  to  describe  the 
vibrational  motion  of  a  diatomic  molecule.  We  consider  a  diatomic  molecule 


with  two  electronic  states,  a  ground  state  |g>  and  an  excited  state  |e>. 
Its  Hamiltonian  is 


H  -  |g>Hg<g|  +  |c>(Wge  +  He)<e|  .  (48) 

We  assume  the  ground  state  potential  to  be  harmonic  and  the  excited  state 
potential  to  be  given  by  a  Morse  oscillator.  We  then  have 


H  ■  §-  +  Jyw2(x-x  )2 
g  Zy  g  g 

2  -g(x-x  )  , 

H  -  +  W  +  D  (l-e  ) 

e  2y  g,e  e 


(49a) 


(49b) 


where  \i  is  the  reduced  mass  of  the  molecule,  u>  is  the  vibrational  frequency 
on  the  lower  potential  surface,  W  is  the  excitation  energy  from  lower  to 
the  upper  surface,  D  is  the  equilibrium  dissociation  energy  of  the  upper 
potential  surface,  and  p  is  a  constant  given  by  [21] 


»  •  'Tip'1 


(49c) 


where  c  is  the  velocity  of  light,  h  is  Planck's  constant,  and  u>  is  the 
vibrational  frequency  that  the  anharmonic  oscillator  would  have  classically 
for  an  infinitesimal  amplitude.  For  our  purposes  we  assume 

w  *  4395  cm  p  *  0.5  a.u.,  D  »  38,310  cm  * ,  D  *  hcB 
e  e  e  e 


P  -  1.93  A  ,  x  -  0.504  A,  W 

g  g»e 


26,230  cm*1,  x  -  0.6325  A  .(50) 


We  consider  the  molecule  initially  to  be  in  its  ground  vibrational  state 
satisfying  the  minimum  uncertainty  condition 


and  we  set  initially  o.  *  0.604  A  and  *  0.0  gm  cm/s.  We  now  assume  that 
at  time  t  *  0  there  is  a  Franck-Condon  transition  from  the  ground  to  the 
excited  potential  surface.  After  this  transition  the  molecule  will  start 
executing  vibrational  motion  about  the  excited  state  equilbrlum  position  X  • 
To  study  this  vibrational  motion,  we  solve  the  TDSCF  set  of  REM  (34)  fn 
dimensionless  form,  where  the  dimensionless  quantities  o  '  s  are  related  to 
the  o  's  as 


/mjL)  °2 9 


mhu>  °4 


where  u>  -  2*cu>  with  initial  conditions  (50),  o-  ■  5.0,  o,  -  0.05  and  with 

up  to  200  time  steps  on  the  order  of  -0.3x10  s7  Variations  of  o^(t)  and 

o^Ct)  with  time  are  shown  in  Figs*  1  and  2,  respectively.  We  have  used  the 
ordinary  differential  equation  solver  technique  of  Gear  [22]  to  solve  the 
TDSCF  set  of  equations  (34).  In  Fig.  3  we  elaborate  further  on  the 

performance  of  our  SCF  approach  by  tracing  the  path  of  o.  over  the  excited 

state  potential  surface.  As  evident  from  Figs.  1  and  i,  given  the  initial 
and  o~  on  the  potential  surface,  which  for  the  present  case  is  o.  -  0.604 
A  and  d2  *  0.0  gm  cm/s,  our  TDSCF  method  describes  enharmonic  vibrational 
motion  of^the  diatomic  molecule  over  this  surface  from  of  ■  0.604  A  to  1.08 


For  the  sake  of  comparison,  we  also  solve  in  dimensionless  form  the 
TDSCF  set  of  REM  (35),  which  describe  the  variations  of  o^CtVs  with  time 
under  the  quadratic  potential  approximation.  We  use  the  same  set  of  initial 
conditions  as  above.  Time  variations  of  o.(t)  and  o-(t)  for  the  present 
case  are  shown  in  Figs.  4  and  5,  respectively.  In  Fig.  J,  we  trace  the  path 
of  o(t)  obtained  under  the  quadratic  potential  approximation  using  a  dashed 
line  to  illustrate  the  performance  of  Hellers  method  compared  to  our  TDSCF 
method. 
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Fig.  1.  The  dimensionless 
mean  displacement  6^  vs 
time  for  a  Gaussian  wave- 
packet  propogated  on  a 
Morse  potential. 

Fig.  2.  The  dimensionless 
mean  momentum  82  vs  time 
for  a  Gaussian  wavepacket 
on  a  Morse  potential. 

Fig.  3.  The  Morse  (solid 
line)  potential  function  of 
an  arbitrary  diatomic 
molecule.  The  dashed  line 
is  obtained  by  considering 
up  to  the  quadratic  terms 
in  o.(t)  in  the  Morse 
function. 
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Fig.  4.  The  dimensionless  mean  dis¬ 
placement  6^  vs  time  for  a  Gaussian 
wavepacket  propogated  on  a  locally 
quadratic  form  of  the  Morse  potential 


Fig.  5.  The  dimensionless  mean  momentum 
vs  time  for  a  Gaussian  wavepacket 
propogated  on  a  locally  quadratic  form 
of  the  Morse  potential. 


VI.  DISCUSSION 


We  have  described  a  method  for  studying  the  dynamical  properties  of 
irreversible  statistical  systems.  Irreversibility  is  introduced  into  our 
system  through  quantum  measurements  [20],  and  this  enables  us  to  make  use  of 
the  maximum  entropy-based  formulation  (MEF).  Use  of  MEF  in  constructing  the 
reduced  density  matrix  (2)  eliminates  the  necessity  of  performing  tedious 
thermal  averaging  [9].  Therefore,  the  present  TDSCF  method  will  be 
particularly  suitable  for  studying  the  various  dynamical  processes  in 
condensed  phases.  The  present  development  resembles  the  derivation  of 
thermodynamic  theorems  from  statistical  mechanics  due  to  the  fact  that  the 
construction  of  the  density  matrix  and  the  corresponding  REM  are  independent 
of  the  specific  nature  of  the  Hamiltonian.  For  this  reason,  we  find  that 
the  present  method  can  be  used  to  describe  the  relaxation  of  a  system  to 
thermal  equilibrium  with  a  thermal  bath  under  the  exact  and  the  quadratic 
potential  approximations.  Under  certain  conditions  (20),  the  present  method 
can  also  be  used  to  describe  the  time  evolution  of  systems  in  pure  states. 
The  derivation  of  the  REM  are  based  on  a  projection  scheme,  and  the 
projection  operators  are  defined  in  terms  of  the  MEF  density  matrix.  The 
TDSCF  set  of  equations  (35)  and  (40),  which  describe  the  time  evolution  of 
pure  states,  have  been  shown  to  be  quite  useful  for  describing  a  variety  of 
molecular  dynamical  processes,  including  molecular  scattering,  electronic 
spectra,  dissociation  of  clusters  and  thermal  desorption  from  surfaces 
[9,23].  The  present  phase  space  TDSCF  method  enjoys  all  these  advantages. 

In  deriving  the  TDSCF  set  of  equations,  we  have  not  had  to  make  the 
assumption  that  the  exact  nonequilibrium  statistical  density  is  in  some 
sense  approximately  equal  to  the  local  equilibrium  one,  and  thus  the  present 
method  is  much  more  general  than  the  local  equilbirium  formulations.  A 
close  look  at  our  TDSCF  set  of  equations  (18)  shows  that  they  do  not  contain 
H.  That  is,  even  though  we  started  our  development  using  the  quantum 
Liouville  equation,  the  time  evolution  of  our  MEF-based  density  functions 
(28)  and  (29)  is  described  by  a  classical  TDSCF  set  of  REM.  Therefore,  the 
present  MEF-based  TDSCF  method  is  completely  classical.  This  in  turn 
suggests  that  the  present  procedure  may  be  repeated  for  classical  mechanics 
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by  replacing  L  in  Eq.  (5)  with  the  classical  Liouville  operator.  Each 
single-particle  density  function  $,(x,#x';t)  should  then  be  replaced  by  a 
phase  space  distribution  which  is  Gaussian  in  x.  and  p,  (28),  We  can  then 
repeat  the  present  procedure  to  obtain  the  TDSCF  set  of  REM  (18),  and  hence 
to  confirm  their  classical  nature.  Our  TDSCF  method  represents  the  lowest 
order  of  a  systematic  expansion,  (14)  and  (15),  and  may  therefore  be 
improved  by  incorporating  correlation  terms  order  by  order.  Inclusion  of 
the  correlation  terms  will  cause  our  REM  to  contain  Hi,  and  hence  will  depart 
from  the  classical  picture.  Therefore,  the  correlation  terms  may  be 
considered  as  quantum  corrections  to  our  classical  description  [19]. 
However,  for  harmonic  systems  with  normal  mode  x  's,  the  TDSCF  set  of  REM 
are  exact.  An  alternative  way  to  improve  our  TDSCF^description  would  be  to 
include  cubic  and  higher  moments  to  construct  each  single-particle  density 
matrix  ip ,  (x «  ,x ' ;  t) .  This  would  then  be  a  departure  from  the  Gaussian 
picture. J  J  J 

Although  the  inclusion  of  the  correlation  terms,  (1A)  and  (15),  and  the 
higher  moments,  (2),  would  improve  our  TDSCF  description,  the  product 
ansatz,  (1),  for  the  N-particle  density  function  implies  neglect  of  exchange 
effects  and  an  incomplete  account  of  quantum  mechanical  correlations.  This 
is  one  of  the  limitations  of  our  single-particle  description  of  an 
interacting  N-particle  statistical  system.  Implementation  of  the  exchange 
effects  for  equilibrium  Bose  and  Fermi  systems  are  available  in  the 
literature  [4,24],  Again,  the  present  development  is  restricted  to  one 
dimensional  phase  space.  Extension  to  the  simulation  of  equilibrium  and 
nonequilibrium  statistical  systems  in  three-dimensional  phase  space  will  be 
reported  in  the  future, 
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